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Abstract 

For the design and optimization of CO2 recovery from alcoholic fermentation pro- 
cesses by distillation, models for vapor-liquid equilibria (VLE) are needed. Two such 
thermodynamic models, the Peng-Robinson equation of state (EOS) and a model 
based on Henry's law constants, are proposed for the ternary mixture N2-|-02-|-C02. 
Pure substance parameters of the Peng-Robinson EOS are taken from the litera- 
ture, whereas the binary parameters of the Van der Waals one-fluid mixing rule 
are adjusted to experimental binary VLE data. The Peng-Robinson EOS describes 
both binary and ternary experimental data well, except at high pressures approach- 
ing the critical region. A molecular model is validated by simulation using binary 
and ternary experimental VLE data. On the basis of this model, the Henry's law 
constants of N2 and O2 in CO2 are predicted by molecular simulation. An easy-to- 
use thermodynamic model, based on those Henry's law constants, is developed to 
reliably describe the VLE in the C02-rich region. 
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1 Introduction 



The recovery of CO2, produced during alcoholic fermentation, is an economically and 
ecologically interesting process for breweries. CO2, which is the main component of the 
fermentation flue g be collected and purified to be used as an auxiliary material 

during the production process or for other purposes. E.g., CO2 is used to avoid contact of 
the beer product with O2 from the atmosphere to minimize oxidation processes and flavor 
derogation. It is needed in this case with a high purity, particularly the O2 concentration 
has to be below 5 ppm [H, [2]. 

CO2 originates from the anaerobic metabolism of yeast during fermentation processes. 
Yeast, mostly saccharomyces cerevisiae uvarum varians carlsbergensis, ferments the hex- 
oses glucose and fructose as well as the disaccharides saccharose and maltose and the 
trisaccharides maltotriose. The fermentation produces per liter of beer about 42 g CO2. 
Wort solves about 4 g/1, thus 38 g/1 are released and may predominantly be recovered. 
Fermentation by-products and components from the atmosphere in the fermentation tank 
contaminate the emerging CO2 and can be separated by different cleaning steps via acti- 
vated carbons and silica gels [3]. The permanent gases N2 and O2 from the atmosphere, 
however, are the main obstacles in the recovery of CO2. Available experimental data in- 
dicates low solubilities of O2 and N2 in liquid CO2 at temperatures below —35 °C, which 
allow a sufficiently high purification [1]. 

N2 and O2 can be separated from CO2 by distillation, recovering up to 66 % (i.e. 25 
g/1) of the produced CO2. Common one-stage cooling devices used in breweries, mostly 
with NH3 as working agent, are unable to provide the required low temperatures down to 
—55 °C [5]. Therefore, an additional low temperature stage has to be considered. CO2 
itself is an appropriate cooling agent for such devices due to its volumetric refrigerating 
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capacity and the optimal pressure range at the required temperatures. In a pilot plant 
at the Flensburger Brauerei, Emil Petersen GmbH & Co. KG (Germany), which is dis- 
cussed in [_6j, a cascade cooling device was added. After hquefaction and super-cooling, 
the recovered C02-rich liquid is stored into an interim tank. Subsequently, the liquid 
mixture is heated by three heating sections to remove the super-cooling and to boil out 
the permanent gases N2 and O2. The gaseous fraction is dissipated to avoid a re-solution. 

For the design and optimization of such recovery and cleaning processes, vapor-liquid 
equilibrium (VLE) data for the ternary mixture O2+N2+CO2 at temperatures between 
—55 and —20 °C are needed. The focus of the present work is the assessment of the 
available experimental VLE data and the development of two thermodynamic models for 
that purpose. Particular attention is given to the C02-rich region. Beside the classi- 
cal approach with the Peng- Robinson equation of state (EOS), molecular modelling and 
simulation were used to develop an easy-to-use model based on Henry's law constants. 
It should be noted that the molecular approach is very much suitable to predict mixing 
properties for a wide variety of fluids. 

2 Experimental Data 

Experimental VLE data of the ternary mixture N2+O2+CO2 are available only from two 
publications, Zenner et al. [7] and Muirbrook et al. |8j. The later [5] deals exclusively with 
the °C isotherm, thus Zenner et al. [7] is the only ternary source within the regarded 
temperature range. In [7J, the VLE has been measured at the temperature —40.3 °C and 
pressures of 5.17, 6.90 and 12.69 MPa with CO2 liquid mole fractions ranging from 0.73 
to 0.92. At —55 °C pressures of 6.90, 10.35, and 13.10 MPa have been investigated with 
CO2 liquid mole fractions ranging from 0.69 to 0.87. These are high pressure VLE so that 
the very C02-rich region is not covered. 
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Regarding the three binary subsystems, the available experimental data base is much 
better and extensive VLE data can be found. For N2+O2, which consists of the two 
lower boiling components of the present ternary mixture, 15 publications are available 
[8-22]. Both components have very low critical temperatures, i.e. —147.05 °C (N2) and 
—118.57 °C (O2), so that the binary VLE lies considerably lower than the temperature 
range of interest as well. 

The VLE of the binary subsystem N2+CO2 has been investigated in 21 publications 
[6,7,23-41]. CO2 has a much higher critical temperature of 30.98 °C and experimental 
binary VLE data can be found in the regarded temperature range in 12 publications 
[25-36]. 

Also for the third subsystem, i.e. O2-I-CO2, sufficient VLE data is available [6,31,42- 
44]. In this case, numerous data points are known in the relevant temperature range from 
-55 to -20 °C m SSI m] as well. 

In the following, a subset of the experimental data was used to adjust and validate 
the thermodynamic models developed in this work, i.e. [T^ \T7\ UT] for N2-I-O2, [T| for 
N2+CO2, [43] for O2+CO2, and [7J for the ternary mixture N2+O2+CO2. For N2+O2 
and N2-I-CO2 those data sets were selected, which contain a larger number of data points. 

3 Molecular Model 

To describe the intermolecular interactions in the ternary mixture, effective state inde- 
pendent pair potentials were used here, which implies that many-body interactions were 
neglected. For this purpose, the two-centre Lennard- Jones plus point quadrupole (2CL JQ) 
pair potential was employed [16]. It is composed of two identical Lennard- Jones sites a 
distance L apart (2CLJ) and an axial point-quadrupole of momentum Q placed in the 
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geometric centre of the molecule. The intermolecular potential writes as 

12 / \ 6" 



a=l b=l 
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+ (1) 



Herein, Vij is the centre-centre distance vector of two molecules i and j, Tab is one of the 
four Lennard- Jones site-site distances; a counts the two sites of molecule i, b counts those 
of molecule j. The vectors oJi and u}j represent the orientations of the two molecules 
i and j. /q is a trigonometrical function depending on these molecular orientations, cf. 
Gray and Gubbins [TF] . The Lennard- Jones parameters a and e represent size and energy, 
respectively. In total, the 2CLJQ model has four model parameters: a, e, L, and Q. These 
parameters have been adjusted to VLE for numerous pure fluids in prior work [46j. Table 
[1] summarizes the parameters of the three pure fluid molecular models considered here. 

On the basis of existing models for pure fluids, molecular modelling of mixtures reduces 
to specifying the interaction between unlike molecules. Following prior work [48j , a mod- 
ified Lorentz-Berthelot combining rule with one adjustable binary interaction parameter 
^ was used for each unlike Lennard- Jones interaction 

(TAB = , (2) 

Cab = ^- V^A-eB- (3) 

Table [2] summarizes the three binary interaction parameters needed for the ternary mix- 
ture model that were taken from [18]. The interaction between the quadrupolar sites 
is treated in a physically straightforward way without the use of binary parameters. It 
was shown in that this ternary model yields an accurate description of the thermo- 
dynamic properties of this mixture. It can readily be used to predict a wide range of 
thermodynamic properties such as Henry's law constants as discussed below. 
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To calculate the VLE on the basis of this molecular model by simulation, the Grand 
Equilibrium method was applied here. The description of this method can be found 
elsewhere |19] and is not repeated. Only the technical details of the present calculations 
are concisely given in the following. Molecular dynamics simulations for the liquid phase 
were performed in the isobaric isothermal (NpT) ensemble, using Anderson's barostat [50j 
and isokinetic velocity scaling [51] . A total of 864 molecules were placed initially in a fee 
lattice configuration in a cubic simulation volume. Depending upon the density of the 
state point, the reduced membrane mass parameter for the barostat Mm/a^ was chosen 
from 10^^ to 10^^, where m is the molecular mass. The intermolecular interactions were 
evaluated explicitly up to a cut-off radius of 5cr and standard long range corrections were 
used, employing angle averaging as proposed by Lustig [52] . 

For the vapor phase, pseudo-grand canonical (p-/i\^T) Monte-Carlo simulations were 
performed. The cut-off radius was also Tc = 5a and the long range corrections were 
considered. The maximum displacement was set to 5% of the simulation box length, which 
was chosen to yield on average 300 to 400 molecules in the volume. After 1 000 initial 
cycles in the canonical (NVT) ensemble starting from a fee lattice, 9 000 equilibration 
cycles in the p-fiVT ensemble were performed. One cycle is defined here to be a number 
of attempts to displace and rotate molecules equal to two times the actual number of 
molecules plus three insertion and three deletion attempts. The length of the production 
run was 100 000 cycles. In this way, VLE data on the basis of the molecular model 
were generated. The results from simulation are presented in section [H] and validated by 
comparison to experimental data. 
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4 Peng-Robinson Equation of State 



Cubic EOS offer a compromise between generality and simplicity that is suitable for 
numerous purposes. They are excellent tools to correlate experimental data and are 
therefore often used for many technical applications. In the present work, the Peng- 
Robinson EOS with the Van der Waals one-fluid mixing rule was adjusted to binary 
experimental data and validated regarding the ternary mixture. The Peng-Robinson 
EOS [Mj is defined by 

RT a 



V — b v{v + h) + h{y — h) ' 
where the temperature dependent parameter a is defined by 



(4) 
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0.45724 



Pc 



1 + (0.37464 + 1.54226 u - 0.26992 cj^) 1 - W — 
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The constant parameter h is 



0.07780 



RTr 



Pc 



(6) 



Therein, critical temperature Tc, critical pressure pc, acentric factor tu, and the ideal gas 
constant R of the pure substances are needed, cf. Table [31 The values were taken from 



To apply the Peng-Robinson EOS to mixtures, mixed parameters am and hm have to 
be defined. For this purpose a variety of mixing rules can be found in the literature. 
Here, the Van der Waals one-fluid mixing rule [22] was chosen. It defines the temperature 
dependent parameter of the mixture as 



(7) 



The indices i and j denote the components, with 



where kij is an adjustable binary parameter. The constant parameter of the mixture is 
defined as 



This classical thermodynamic model was used to fit the experimental VLE data of the 
three binary subsystems, i.e. N2+O2, O2+CO2, and N2+CO2. The three binary param- 
eters kij were adjusted to the same experimental data as the binary parameters of the 
molecular model ^. It turned out that temperature independent kij values are sufficient 
in all cases in the regarded range of states, cf. Table |H The results of the Peng- Robinson 
EOS are presented in section [6] and validated by comparison to experimental data. 

5 Henry Model 

Molecular simulation allows the prediction of Henry's law constants on the basis of a 
given molecular model straightforwardly. Different approaches have been proposed in the 
literature, e.g. [551 [56] . The Henry's law constant is related to the residual chemical 
potential of the solute i at infinite dilution [57] 



where is the Boltzmann constant, T the temperature, and p the density of the solvent. 

For the calculation of Henry's law constants, a series of simulations, ranging from to 
—55 °C, were performed with at intervals of 5 °C. To evaluate Pi°°, Widom's test insertion 
method [58] was used. Therefore, 3456 test molecules representing the solute i were 
inserted after each time step at random positions into the liquid solvent and the potential 
energy between the solute test molecule and all solvent molecules ipi was evaluated within 
the cut-off radius 
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where the brackets represent the ensemble average in the NpT ensemble. Note that 
appropriate long range corrections [51] have to be applied. The residual chemical potential 
at infinite dilution from this procedure is solely attributed to the unlike solvent-solute 
interactions. The mole fraction of the solute in the solvent is exactly zero, as required for 
infinite dilution, since the test molecules are instantly removed after the potential energy 
calculation. Simulations were performed at a specified temperature and the according 
pure substance vapor pressure of CO2. The results for the Henry's law constants Hi are 
given as functions of the temperature as shown in Figure [TJ Linear functions were found 
to be sufficient to fit the data. The resulting equations are for N2 in CO2 

i/N2/MPa = 178.21 - 0.47998 T/K, (12) 

and for O2 in CO2 

ii'o2/MPa = 104.26 - 0.23214 T/K. (13) 

It can be seen in Figure [1] that the predicted Henry's law constant of O2 in CO2 agrees 
well with the experimental data [HIMIIIS], especially in the low temperature region. But 
a considerable scatter of experimental data has to be noted. For N2 in CO2 systematic 
deviations between simulation and experiment were found. At —55 °C the agreement is 
very good, but with increasing temperature the data sets diverge, where simulation yields 
higher values. The deviation is 14% at —40.3 °C. 

The classical approach to model VLE on the basis of Henry's law constants includes 
the activity coefficient 7^ and the fugacity coefficient 0j. The phase equilibrium condition 
for the two low boiling components i = N2, O2 is then given by [53] 

Hi exp J f v°°dp \ Xi-fi=pyi (j)i. (14) 

I JPC02 J 

Here, Xj and i/i are the mole fractions in the saturated liquid and vapor, respectively, and 
is the partial molar volume of the solute at infinite dilution. The exponential term, 
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known as Krichevski-Kasarnoski correction, accounts for the higher pressure of the hquid 
mixture compared to the pure solvent vapor pressure. In the pressure range of interest, 
its influence is small so that this correction term was set to unity. 

For the solvent CO2, the equilibrium condition includes the pure substance vapor 
pressure p*co2 instead of Henry's law constant 

P'C02 XC02 7C02 = P I/CO2 0CO2- (15) 

Therefore, a correlation for p'^co2 [59] was taken from the literature 

4 

Hp'co2/Pc) = Y.M^ - {T/Tjr/{T/T,), (16) 

i=l 

where the parameters Ai and rii are listed in Table O 

With the previously adjusted Peng-Robinson EOS, activity coefficients were calculated 
in the temperature and composition range of interest. It was found that they are between 
1 and 1.001 and were thus set to unity in the present Henry model. 

The fugacity coefficients, which describe the non-ideality of the vapor phase, were 
calculated with the Peng-Robinson EOS as well. The values do have a considerable 
temperature dependence, cf. Figure [2], and were thus included in the present Henry 
model. Quadratic fits were found to be sufficient for all components, i.e. for N2 

0N2 = 5.5276 - 0.038326 T/K + 0.00008294 (T/K)^ (17) 

for O2 

^02 = 0.54440 - 0.00053697 T/K + 0.00001082 (r/K)^, (18) 

and for CO2 

0CO2 = 0.51216 + 0.0046540 T/K - 0.00001082 {T/Kf. (19) 

Equations (12) to (19) define the present Henry model to describe the VLE of the ternary 
mixture in the C02-rich region. 
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6 Results and Discussion 



In this section, present thermodynamic models are vahdated against experimental data. It 
is started with the three binary subsystems and subsequently the ternary case is regarded. 
Particular attention is given to the molecular model as it was used to predict data for 
development and validation of the Henry model. 

6.1 Nitrogen+ Oxygen 

In Figure [31 molecular simulation results and the Peng-Robinson EOS are compared 
to experimental VLE data [121 El EI] of N2+O2 at —153.15 °C. This temperature is 
considerably lower than the target range of —55 to —20 °C due to the fact that both 
components have very low critical temperatures. It can be seen that the Peng-Robinson 
EOS correlates the experimental data well. By simulation, an equimolar composition 
in the liquid phase was regarded, where the mixing effect is strongest. The agreement 
between the three data sets is very good, which is also the case for other temperatures 
(not shown here). 

6.2 Nitrogen+ Carbon Dioxide 

Figure m depicts the VLE of N2+CO2 at —40.3 °C including experimental data [7], molec- 
ular simulation results, Peng-Robinson EOS, and Henry model. The experimental data 
shows some scatter, the simulation results are within this error bound. This also holds 
for the Peng-Robinson EOS for pressures up to 8 MPa. Approaching the critical region 
of the mixture, it is found that the Peng-Robinson EOS overshoots considerably which 
is a well known problem of this thermodynamic model. The Henry model is insufficient 
at °C (not shown here), but for lower temperatures, e.g. cf. Figure HI it agrees well in 
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the C02-rich region. By closer inspection of the data, which will be made below, some 
deviations are found on the dew line that can hardly be seen with the resolution chosen 
for Figure |H 

6.3 Oxygen+ Carbon Dioxide 

Figure [5] presents the VLE data of O2+CO2 at —40.3 °C from experiment [13], molecular 
simulation, Peng-Robinson EOS, and Henry model. As before, the experimental data 
shows some scatter, particularly on the dew line. Again, it can be seen that the Peng- 
Robinson EOS overshoots in the critical region of the mixture. The agreement between 
Peng-Robinson EOS and Henry model is good in the C02-rich region. The molecular 
model shows reliable results, also for other temperatures (not shown here). As for the 
previous binary mixture, some deviations are observed on the dew line. 

6.4 Nitrogen+Oxygen+ Carbon Dioxide 

In Figures [6] and [7| simulation results and the Peng-Robinson EOS are compared with 
experimental ternary VLE data of the ternary system at —40.3 and —55 °C [7J. Due to 
the fact that experimental data is available only at high pressures between 5.17 and 13.10 
MPa, the Henry model is not included into this validation. The Peng-Robinson EOS 
represents the VLE well at pressures of 6.9 MPa and below as depicted in Figures [6] and 
[71 But it shows larger deviations at higher pressures approaching the critical region (not 
shown here). In the target region of state points, molecular simulation data is throughout 
in very good agreement with the experimental data. 
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6.5 Carbon Dioxide-Rich Region 

The present Henry model is developed for the design of technical applications in the 
C02-rich composition range. Therefore, further validations in this region, at state points 
where no experimental data is available, are presented here. The molecular model, being 
validated on the basis of experimental mixture VLE as discussed above, was used as the 
benchmark here. 

Figure [8] presents the pressure over the vapor mole fractions of N2 and O2 in VLE 
for constant liquid mole fractions XN2=a^O2=0.01 in the temperature range from —55 to 
—20 °C. Simulation data, Peng-Robinson EOS, and Henry model are compared. It can be 
seen that simulation data and Peng-Robinson EOS agree very well, deviations are minor. 
The Henry model deviates somewhat, yielding approximately 8% too low pressures and 
5% too low vapor mole fractions. 

The limits, where the Henry model shows deviations of less than 2% from molecular 
model and Peng-Robinson EOS, were investigated. Table [6] shows that the Henry model 
is reliable for CO2 liquid mole fractions above 0.995. This limit is examined in Figure [H 
where the pressure over the vapor mole fractions of N2 and O2 in VLE for constant liquid 
mole fractions XN2=a;o2=0.0025 is shown in the temperature range from —55 to —20 °C. 
The results from all three models agree well, proving the reliability of the thermodynamic 
models presented in this work. 

7 Conclusion 

In this work three thermodynamic approaches, i.e. a molecular model, the Peng-Robinson 
EOS, and a model based on Henry's law constants, were used to investigate the VLE 
of the ternary mixture N2+O2+CO2 in the temperature range from —55 to —20 °C. A 
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thorough vahdation by comparison to experimental VLE data was made, where possible. 
The molecular model and the Peng-Robinson EOS are appropriate throughout, except 
in the critical region. For the very C02-rich region, which is important for purification 
processes, the computationally convenient Henry model can be used reliably. 
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List of symbols 



Latin Letters 





parameter of Peng-Robinson equation of state 


A 


coefficients of correlation 


b 


parameter of Peng- Robinson equation of state 


fq 


trigonometrical function depending on molecular orientations 


H 


Henry's law constant 


i 


molecule index 


j 


molecule index 


ks 


Boltzmann's constant, ks = 1.38066-10^^ J/K 


kij 


binary parameter of Peng-Robinson equation of state 


L 


elongation 


Tfl 


molecular mass 


M 


membrane mass parameter 


ft 


exponent of correlation 


P 


pressure 


Q 


quadrupolar momentum 


R 


ideal gas constant 


T 


temperature 


u 


pair potential 


V 


molar volume 


V 


extensive volume 


X 


mole fraction in liquid phase 


y 


mole fraction in vapor phase 
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Greek Letters 



7 activity coefficient 

e Lennard- Jones energy parameter 

^ binary interaction parameter 

II chemical potential 

p number density 

a Lennard- Jones size parameter 

(f) fugacity coefficient 

ip potential energy of test molecule 

cu acentric factor 



Subscripts 



a count variable for molecule sites 

A related to component A 

b count variable for molecule sites 

B related to component B 

c critical value 

i related to component i 

ij related to components i and j 

j related to component j 

m mixture 

Q quadrupole 
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Superscripts 

oo at infinite dilution 

s saturated 

Abbreviations 



2CLJ two-center Lennard- Jones 

2CLJQ two-center Lennard- Jones plus point quadrupole 

EOS equation of state 

NpT isobaric-isothermal ensemble 

NVT canonic ensemble 

IjlVT grand canonical ensemble 

VLE vapor-liquid equilibria 



Vector properties 

Vij center-center distance vector between two molecules i and j 

uj orientation vector of a molecule 
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Table 1: Parameters of the molecular models for the pure fluids, taken from [46j. 



Fluid 


a/k 


(e/^B) /K 


L/k 


g/DA 


N2 


3.3211 


34.897 


1.0464 


1.4397 


O2 


3.1062 


43.183 


0.9699 


0.8081 


CO2 


2.9847 


133.22 


2.4176 


3.7938 



Table 2: Binary interaction parameters of the molecular model, taken from [18]. 



Mixture 




N2 + O2 


1.007 


N2 + CO2 


1.041 


O2 + CO2 


0.979 



Table 3: Pure substance parameters of the Peng- Robinson EOS, taken from [54j . 



Fluid 




Pc/MPa 


UJ 


N2 


126.19 


3.3958 


0.0372 


O2 


154.58 


5.0430 


0.0222 


CO2 


304.13 


7.3773 


0.2239 
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Table 4: Binary parameters of the Van der Waals one- fluid mixing rule adjusted in the 
present work. 



Mixture 


k- ■ 

1-3 


N2 + O2 


-0.0119 


N2 + CO2 


0.0015 


O2 + CO2 


0.124 



Table 5: Parameters of the vapor pressure correlation for CO2, taken from [59] 



i 


A, 


Hi 


1 


-6.95626 


1 


2 


1.19695 


3/2 


3 


-3.12614 


3 


4 


2.99448 


6 



Table 6: Minimum mole fractions of CO2 in the liquid and vapor where the Henry model 
deviates by less than 2 % from the molecular model and the Peng-Robinson EOS. 



T/°C 


XC02 


2/C02 


-20 


0.994 


0.872 


-35 


0.995 


0.811 


-55 


0.996 


0.706 
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